## This code creates Table A6 and Figure A1.

neighbor_dist_cat <- seq(10,100,10) # Town distance cutoff points

data <- colonial_towns_shp[!duplicated(colonial_towns_shp$latitude),]

data <- colonial_towns_shp[colonial_towns_shp$year == 1912,]

nb<- tri2nb(cbind(data$latitude,data$longitude))

lw <- nb2listw(nb)

# Results for Table A6

moran.test(data$coercive,lw, alternative="greater")

moran.test(data$m_police,lw, alternative="greater")

# Figure A1 (a)

coercive_spatial_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  data <- colonial_towns_shp[!duplicated(colonial_towns_shp$latitude) & colonial_towns_shp$neighbor_dist<=neighbor_dist_cat[j],]
  nb<- tri2nb(cbind(data$latitude,data$longitude))
  lw <- nb2listw(nb)
  model <- spatialreg::lagsarlm(coercive~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data,lw, tol.solve=1.0e-30,zero.policy = TRUE)
  coercive_spatial_BAGO[j,"coefficient"] <- model$coefficients[[2]]
  coercive_spatial_BAGO[j,"se"] <- summary(model)$Coef[,2][[2]]
  coercive_spatial_BAGO[j,"upper"] <- coercive_spatial_BAGO[j,"coefficient"] + (1.96*coercive_spatial_BAGO[j,"se"])
  coercive_spatial_BAGO[j,"lower"] <- coercive_spatial_BAGO[j,"coefficient"] - (1.96*coercive_spatial_BAGO[j,"se"])
  coercive_spatial_BAGO[j,"upper_10pc"] <- coercive_spatial_BAGO[j,"coefficient"] + (1.645*coercive_spatial_BAGO[j,"se"])
  coercive_spatial_BAGO[j,"lower_10pc"] <- coercive_spatial_BAGO[j,"coefficient"] - (1.645*coercive_spatial_BAGO[j,"se"])
  coercive_spatial_BAGO[j,"pch"] <- ifelse(sign(coercive_spatial_BAGO[j,"upper"])==sign(coercive_spatial_BAGO[j,"lower"]),"black",ifelse(sign(coercive_spatial_BAGO[j,"upper_10pc"])==sign(coercive_spatial_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/coercive_spatial_BAGO.png",sep=""), width=6.5, height=6, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     main="Number of civilian and military police stations (count)",
     xlab="Distance from a historical location (km)",
     ylab="Effect of early centralization")
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,coercive_spatial_BAGO$lower/2,x.axis,coercive_spatial_BAGO$upper/2)
points(x.axis,coercive_spatial_BAGO$coefficient/2, pch=21, bg=coercive_spatial_BAGO$pch)
abline(h=0, lty=3)
dev.off()

# Figure A1 (b)

m_police_spatial_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  data <- colonial_towns_shp[!duplicated(colonial_towns_shp$latitude) & colonial_towns_shp$neighbor_dist<=neighbor_dist_cat[j],]
  nb<- tri2nb(cbind(data$latitude,data$longitude))
  lw <- nb2listw(nb)
  model <- spatialreg::lagsarlm(m_police~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data,lw, tol.solve=1.0e-30,zero.policy = TRUE)
  m_police_spatial_BAGO[j,"coefficient"] <- model$coefficients[[2]]
  m_police_spatial_BAGO[j,"se"] <- summary(model)$Coef[,2][[2]]
  m_police_spatial_BAGO[j,"upper"] <- m_police_spatial_BAGO[j,"coefficient"] + (1.96*m_police_spatial_BAGO[j,"se"])
  m_police_spatial_BAGO[j,"lower"] <- m_police_spatial_BAGO[j,"coefficient"] - (1.96*m_police_spatial_BAGO[j,"se"])
  m_police_spatial_BAGO[j,"upper_10pc"] <- m_police_spatial_BAGO[j,"coefficient"] + (1.645*m_police_spatial_BAGO[j,"se"])
  m_police_spatial_BAGO[j,"lower_10pc"] <- m_police_spatial_BAGO[j,"coefficient"] - (1.645*m_police_spatial_BAGO[j,"se"])
  m_police_spatial_BAGO[j,"pch"] <- ifelse(sign(m_police_spatial_BAGO[j,"upper"])==sign(m_police_spatial_BAGO[j,"lower"]),"black",ifelse(sign(m_police_spatial_BAGO[j,"upper_10pc"])==sign(m_police_spatial_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/m_police_spatial_BAGO.png",sep=""), width=6.5, height=6, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     main="Number of military police stations (count)",
     xlab="Distance from a historical location (km)",
     ylab="Effect of early centralization")
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,m_police_spatial_BAGO$lower/2,x.axis,m_police_spatial_BAGO$upper/2)
points(x.axis,m_police_spatial_BAGO$coefficient/2, pch=21, bg=m_police_spatial_BAGO$pch)
abline(h=0, lty=3)
dev.off()


